#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include<time.h>
#include<mpi.h>

/*
*   Nesta implementação, utilizamos a matriz de coeficientes e 4 vetores,
*   X que guarda os valores dos coeficientes na iteração k+1
*   Xant que guarda os valores dos coeficientes na iteração k
*   B que guarda os valores de B*(B dividido pelo elemento da diagonal principal)
*   row_teste que armazena os coeficientes da linha a ser testada no final
*
*   row_teste foi necessario pois alteramos a matriz incial para efetuar os calculos
*
*   optamos por nao guardar o resultado de todas as iterações para otimizar o uso
*   de memória, tendo em vista matrizes grandes.
*
*/

float *aloca_vetor(int ordem);
void preenche_row_teste(int ordem,int row,float **matriz,float *row_test);

int malloc2dfloat(float ***array, int n, int m) {

    /* allocate the n*m contiguous items */
    float *p = (float *)calloc(n*m,sizeof(float));
    //float *p = (float *)malloc(n*m*sizeof(float));
    if (!p) return -1;

    /* allocate the row pointers into the memory */
    (*array) = (float **)malloc(n*sizeof(float*));
    if (!(*array)) {
        free(p);
        return -1;
    }

    /* set up the pointers into the contiguous memory */
    int i;
    for (i=0; i<n; i++)
        (*array)[i] = &(p[i*m]);

    return 0;
}

int free2dfloat(float ***array) {
    /* free the memory - the first element of the array is at the start */
    free(&((*array)[0][0]));

    /* free the pointers into the memory */
    free(*array);

    return 0;
}





/*
*   Funcao que efetua a leitura dos valores da matriz de coeficientes no arquivo
*/
int preenche_matriz(FILE *fp,int ordem,float **matriz){
    int i,j;
    float valor;
    for(i=0;i<ordem;i++){
        for(j=0;j<ordem;j++){
            fscanf(fp,"%f",&valor);
            matriz[i][j] = valor;
        }
    }
    return 0;
}

/*
*   Preenche um vetor para efetuar o teste final
*/
void preenche_row_teste(int ordem,int row,float **matriz,float *row_test){
    int i;
    for(i=0;i<ordem;i++){
        row_test[i]=matriz[row][i];
    }

}

/*
*   Funcao que efetua a leitura dos valores do vetor X no arquivo
*/
int preenche_B(FILE *fp,int ordem,float *B){
    int i;
    float valor;
    for(i=0;i<ordem;i++){
        fscanf(fp,"%f",&valor);
        B[i] = valor;
    }
    return 0;
}

/*
*   Faz o calculo de matriz*, X0 e B*
*/
int calc_mX(int ordem,float **matriz,float *X,float *B, float *Xant){
    int i,j;
    float div;
    //loop para dividir os membros da linha pelo valor da diagonal principal da linha
    for(i=0;i<ordem;i++){
        div = matriz[i][i];
        for(j=0;j<ordem;j++){
            matriz[i][j] /= div;
        }
        Xant[i] = B[i]/div;
    }
    return 0;
}

//void prepareMatrices(float **mA, float Xant, int ordem, int ordem){
void prepareMatrices(float **mA, float *mB, int lines, int colums){
    int i = 0, j = 0;
    for(i=0;i<lines;i++){ //sao excluidos os elementos das diagonales
        for(j=0;j<colums;j++){
            if(i!=j){
                mA[i][j] = -(mA[i][j]); //pela forma da equacao os elementos mudao de signo
            }
            else{
                mA[i][j] = mB[i]; //elemento do vetor b
            }
        }
    }
}

/*
*   Funcao que calcula o valor de Xk+1 e atualiza o valor de Xk
*/
void processingMatrix(int base, int limite, int ordem, float **main, float *X, float *Xant ){
    int i, j;
    for (i = base; i < limite; i++) {

        X[i] = 0;
        for (j = 0; j < ordem; j++) {
            if(i!=j){
                X[i] += main[i][j] * Xant[j];
            }else{
                X[i] += main[i][j];
            }
        }
    }
}

/*
*   Retorna o erro entre as iteracoes
*/
float relativeError(int lines, float *X, float *Xant){
    int i;
    float maxD = -3000;
    float maxN = -3000;
    float tmpd = -3000;
    float tmpn = -3000;
    for (i = 0; i < lines; i++) {
        tmpn = fabs(X[i] - Xant[i]);
        tmpd = fabs(X[i]);
        if(tmpn > maxN)
            maxN = tmpn;
        if(tmpd > maxD)
            maxD = tmpd;
    }

    return maxN/maxD;
}



/*
*   Main
*/

int main(int argc, char *argv[]){

    FILE *fp;
    int ordem, fila, iteracoes;//dados do arquivo
    int i, j, k, h;//iterador
    int base, limite; //bloco de posicoes a ser computado
    float tol,erro = 99,valorB;//tolerancia e erro,valorB armazena o resultado da equacao original
    float **matriz;//matriz
    float *X, *B, *Xant, *row_test, *Xaux;//vetores
    float result=0;
    clock_t inicio, fim;
    double tempo;
    // variaveis usadas para MPI
    int numprocs, rank, namelen, rc;
    char processor_name[MPI_MAX_PROCESSOR_NAME];

    if(!(fp = fopen(argv[1],"r"))){//abre o arquivo
        printf("Erro! Verifique se o arquivo existe!\n");
        exit(1);
    }

    //pega dados do arquivo
    fscanf(fp,"%d",&ordem);
    fscanf(fp,"%d",&fila);
    fscanf(fp,"%f",&tol);
    fscanf(fp,"%d",&iteracoes);

    //aloca matriz
    malloc2dfloat(&matriz, ordem, ordem);

    //aloca vetor para efetuar o teste final
    row_test = (float *)malloc(ordem* sizeof(float));

    //aloca vetor de iteracoes
    X = (float *)malloc(ordem* sizeof(float));

    //aloca B
    B = (float*)malloc(ordem* sizeof(float));

    //aloca vetor Xant usado para calcular o erro
    Xant =(float*)malloc(ordem* sizeof(float));

    //aloca vetor Xant usado para calcular o erro
    Xaux =(float*)malloc(ordem* sizeof(float));

    //preenche a matriz com os dados
    preenche_matriz(fp,ordem,matriz);

    preenche_row_teste(ordem,fila,matriz,row_test);

    //preenche o vetor B
    preenche_B(fp,ordem,B);
    valorB = B[fila];

    //calcula matriz*,X0 e B*
    calc_mX(ordem,matriz,X,B, Xant);

    //fecha o arquivo
    fclose(fp);

    //loop de calculo
    inicio = clock();

    //void prepareMatrices(float **mA, float *mB, int lines, int colums){
    prepareMatrices(matriz, Xant, ordem, ordem);

    MPI_Status Stat;
    //inicializa MPI
    rc = MPI_Init(&argc, &argv);

    //verifica se a inicializacao ocorreu
    if(rc != MPI_SUCCESS){
            printf("Falha ao iniciar MPI\n\n");
            exit(-1);
    }
    MPI_Comm_size(MPI_COMM_WORLD, &numprocs);
    MPI_Comm_rank(MPI_COMM_WORLD, &rank);
    MPI_Get_processor_name(processor_name, &namelen);

    for(i=0; i <= iteracoes; i++){
        //se nao for o mestre, faz os calculos para pegar o trecho do vetor X a ser calculado
        if(rank != 0 ){
            base = (rank-1)*(ordem/(numprocs-1));
            limite = rank*(ordem/(numprocs-1));

            //se a ordem nao for divisivel pelo numero de nós
            if(ordem%(numprocs-1) != 0){
                //o ultimo nó calcula as posicoes de X que sobraram
                if(rank == numprocs-1){
                limite = limite + (ordem - limite);
                }
            }
            //efetua o calculo de Xk+1
            processingMatrix(base, limite, ordem, matriz, X, Xant);
            	 
            //envia os dados para o mestre
             MPI_Send(X, ordem, MPI_FLOAT, 0, 1, MPI_COMM_WORLD);
            //recebe dados do mestre para novos calculos e verificar condicao de parada
            MPI_Recv(&erro, 1, MPI_FLOAT, 0, 1, MPI_COMM_WORLD, &Stat);
            MPI_Recv(Xant, ordem, MPI_FLOAT, 0, 1, MPI_COMM_WORLD, &Stat);
            MPI_Recv(X, ordem, MPI_FLOAT, 0, 1, MPI_COMM_WORLD, &Stat);
        }
    		else{//mestre
            //atuzaliza Xant
            for(h = 0; h < ordem; h++){
            Xant[h] = X[h];
            }
            for(k=1; k < numprocs; k++){
            //recebe vetor X dos escravos
            MPI_Recv(Xaux, ordem, MPI_FLOAT, k, 1, MPI_COMM_WORLD, &Stat);
            //faz o calculo de base e limite para saber que posicoes concatenar
            base =(k-1)*(ordem/(numprocs-1));
            limite = (k)*(ordem/(numprocs-1));
            //faz o tratamento para do ultimo escravo calcular as posicoes que sobram
            if(ordem%(numprocs-1) != 0){
                if(k == numprocs-1){
                limite = limite + (ordem - limite);
                }
            }
                //concatena o vetor X pegando as posicoes que interessam do vetor Xaux
                for(h = base; h < limite;h++){
                X[h] = Xaux[h];
                }
            }
            //calcula o erro relativo entre X e Xant
            erro = relativeError(ordem, X, Xant);
            //se o erro for maior que a tolerancia
            if(erro > tol){
                for(j = 0;j < ordem;j++)
                Xant[j] = X[j];
                //envia para todos os escravos o erro, Xant e X para efetuar outra iteracao
                for(k=1;k<numprocs;k++){
                    MPI_Send(&erro, 1, MPI_FLOAT, k, 1, MPI_COMM_WORLD);
                    MPI_Send(Xant, ordem, MPI_FLOAT, k, 1, MPI_COMM_WORLD);
                    MPI_Send(X, ordem, MPI_FLOAT, k, 1, MPI_COMM_WORLD);
                }
            }else{
                break;
            }
        }
    }
    fim = clock();

    result=0;

     if(!rank){ //nó mestre
    printf("Numero de iteracoes %d.\n",(i-1));
    //calcula a aproximacao
        for(i=0;i<ordem;i++){
        result = result + X[i]*row_test[i];
        }
    //calcula tempo utilizado nas iteracoes
    tempo = ((fim-inicio)/(CLOCKS_PER_SEC/1000));
    //imprime resultados
    printf("Tempo utilzado em ms %lf.\n",tempo);
    printf("RowTest %d: => [%f] =? %f\n",fila,result,valorB);
    }
    //mata todos
    MPI_Abort(MPI_COMM_WORLD, MPI_SUCCESS);
    MPI_Barrier(MPI_COMM_WORLD);
    MPI_Finalize();
 
    //libera memoria
    free(matriz);
    free(X);
    free(Xant);
    free(B);

    return 0;
}

